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Abstract 

We measure the critical exponents of the three dimensional Gross-Neveu model with 
two four-component fermions. The exponents are inferred from the scaling behaviour 
of observables on lattice sizes 8 3 , 12 3 , 16 3 , 24 3 , and 32 3 . We find that the model has a 
second order phase transition with v = 1.00(4) and 2 — rj = j/u = 1.246(8). We also 
calculate these exponents, through a second order e-expansion around four dimensions, 
for the three dimensional Higgs- Yukawa model, which is expected to be in the same 
universality class, and obtain 7/1/ = 1.237 and v = 0.948, while recent second order 
l/iV/-expansion calculations give 7/1/ = 1.256 and v = 0.903. We conclude that the 
equivalence of the two models remains valid in 3 dimensions at fixed small Nf values. 
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1 Introduction 



The Gross-Neveu (GN) model describes fermions with four-fermion interaction [|T], H. 
It has a global discrete chiral symmetry, which can break down spontaneously to form 
a chiral condensate. This can be seen as a composite scalar particle that gives a non- 
zero mass for the fermions. Due to its simplicity, the GN model has been studied 
extensively. In two dimensions it is perturbatively renormalisable and asymptotically 
free. In addition, the chiral symmetry is broken. 

In three dimensions (3d) it is renormalisable if the number of flavours, Nf, is large 
enough. Hence, it is also the first renormalisable model known not to be perturbatively 
renormalisable. It has been proven to be renormalisable and non-trivial in dimensions 
between 2 and 4 by means of a l/iVy-expansion B. 

It has been suggested that a model with four-fermion interactions, which leads to 
a spontaneously broken global chiral symmetry with a chiral condensate, could be a 
candidate for the Higgs particle of the standard model 0. Near four dimensions (4d) 
it has been related to Higgs- Yukawa (HY) type models f|, [|. In fact, in 4d both the 
standard electroweak and the GN model (with certain modifications) are trivial and 
can be mapped onto each other J3J and in dimensions 2 < d < 4 the GN model and 
the HY model are equivalent in the framework of l/iVy-expansion [||]. 

The purpose of this work is to enlighten the connection between the composite and 
fundamental Higgs scenarios in a case where the models are not trivial, namely in 
3d. The 1/Nf expansions for the HY model and the GN model are order by order 
equivalent [[J. In order to go beyond the 1/Nf perturbation theory, we ask whether 
this equivalence persists at small Afy's. To answer this we study the critical properties 
of GN model at small Nf by means of a Monte Carlo (MC) simulation. To preserve 
discrete chiral symmetry, the restoration of which we are interested in, we use staggered 
fermions on the lattice and choose the smallest possible value of Nf, that is Nf = 2. A 
finite-size analysis of the numerical simulation gives the critical coupling and critical 
exponent values to be compared to the exponents of the Higgs- Yukawa model. In order 
to obtain a meaningful comparison, one has to go beyond the order e results in the HY 
model . We have therefore extended the calculation and obtain the critical exponents 
to order e 2 . The relevance of the comparison can furthermore be evaluated by using 
recent results from a l/iV/-expansion of the GN model || [7|, ||, |9], [L0| to order 1/Nj. 

The paper is divided into two major parts: section 2 is devoted to the analytical 
results of the Gross-Neveu and Higgs- Yukawa models with the e-expansion while in 
section 3 we describe the Monte Carlo runs performed. Section 4 is devoted to the 
comparison of numerical and analytical results. Within the uncertainties associated on 
the one hand with the statistics of the numerical simulations and on the other hand 
with the still short series expansion in e and 1/Nf, we confirm the equivalence of the 
two 3d models at small Nf. 
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2 Analytical Results 



Here we recall a few known properties of the models under study, and present our new 
analytical calculations, in particular the fixed Nf e 2 expansion of the HY model expo- 
nents. The 1/Nf expansions of the latter are then compared with the 1/Nf expansions 
of the GN model. 

2.1 Continuum Gross-Neveu model 

The continuum GN model with Nf fermion flavours is defined by the Lagrangian 



Nf _ 2 
a=l A 



N f _ 



a=l 

Usually an auxiliary scalar field a is introduced 

Nf 1 



(1) 



£ = £ v> Q w - 9<x)\r{x) - ^ 2 (x), (2) 

a=l Z 



which is formally equivalent to flU) upon integration over a field. 
The GN model has a discrete chiral symmetry 



■ip -> 75^, ip -> —075, cr -> -a, (3) 



which in 3d is spontaneously broken at small couplings. 
The critical exponents for d = 3 to order 1/iVj are || |^, |T 



1 32 64(27tt 2 + 632) 

^~ ~3^iV + 27vr 4 iV 2 ' U 



(5) 



7 _ 64 64(27tt 2 - 304) 
^~ + 3n 2 N + 27tt 4 N 2 ' 
where N = Nf Tr 1 is the total number of fermionic variables. 

2.2 The discretised Gross-Neveu model 

We consider the GN model defined on a 3d symmetric lattice. The discretisation of the 
continuum Lagrangian has to be implemented in such a way as to reproduce the correct 
symmetries in the continuum limit [fTI]. We use staggered fermions, which means that 



in 3d there are 8 doublers. Using 4 component spinors, we can assign the components 
to the corners of the cube, leaving 8/4 = 2 continuum flavours from 1 staggered lattice 
fermion per site. The discretised action reads 

a 2 Nl 

S = N L £^f + E£ Xn (D„ m + £ nm );C (6) 

n ^ n,ma=l 
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The staggered fermion matrix D is given by 

Z j 

where the sum j is over the directions (j = 1,2,3) and rj n j is the staggered fermion 
phase factor 

Vn , = (-l)m+-+n j - K (g) 

with periodic boundary condition in d-1 dimensions and antiperiodic in the last one 

for a finite lattice. 

The mass matrix E is diagonal 

^nm ^nm (9) 

and depends on a field. We choose a discretisation in which a is the average of o at 
the six nearest neighbours of the lattice site n. 

The coupling A is connected to the continuum coupling g by A = g 2 Nf, with Nf = 2Nl 
as explained. 

One may integrate over the Grassmann variables Xx^Xx to express the partition 
function in terms of an effective action 

Se ff = N L E^f-Tr mp + X)]. (10) 

The Nl = oo critical value A°, where chiral symmetry is restored, can be obtained 
from the saddle point equation at a = as 

\° = 0.989. (11) 



Taking into account the quadratic fluctuations (T2| we can obtain the one loop value 
of A c including 1/Nl corrections. As we are interested to rather small Ni values, we 
choose to solve the gap equation of the one-loop effective potential. Solving this gap 
equation on the lattice with linear sizes L = 32 and 64 we estimate 

\l(N L = 1) = 0.800 (12) 

\l(N L = 6) = 0.938 (13) 

This last value agrees with a direct calculation of the perturbative correction in the 
infinite lattice limit. 

2.3 The Higgs- Yukawa model 

The Lagrangian £' of the HY model has a fourth order interaction term and a kinetic 
term added for the a field, 

£'=E^ + 9*W + \™ 2 ° 2 + ~( V) 2 + ^ (14) 

a=l 
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This model becomes renormalisable in four dimensions. It has been argued that the 
terms added are irrelevant for the number of dimensions less than 4 and marginal at 
d=4. Hence, the critical behaviour of HY model and GN model should be identical. 

Due to its renormalisability, the HY model can be studied by means of an e-expansion. 
The two-loop (3 functions and anomalous dimension 7]^ can be obtained from and 
we have computed the mass anomalous dimension rj m . 

Px 4 = -eA 4 + 7 ^(3A2+2^A^ 2 -12^ 4 )+^(^Al-3ArA^ 2 +7iVA 4 ( ? 4 +96iV^), 

(15) 

& = ~\9 + j^s (N/2 + % 3 + _L_(-^±i^^ - 2V + ^hl (16) 
^ = ^ 2 + ^(^-^), (17) 

Vm = -7^p A 4 + (^)l( A 4 + X ± N 9 2 - 2Ng A ) - Va , (18) 

with d = 4 - e and iV = N f Tr 1. 

The corresponding fix points to order e 2 are 



g* 2 1 (N + 66)VN 2 + 132N + 36 - iV 2 + 516iV + 882 



(An) 2 (N + 6) 108(iV + 6) 3 



(19) 



XI _ -N + 6 + VN 2 + 132iV + 36^ 



(4tt) 2 6(iV + 6) 



+ [ -VN 2 + 132iV + 36 (3iV 3 - 43iV 2 - 1545iV - 1224) 

L 2 . (20) 



+37V 4 + 155iV 3 + 2745iV 2 - 2538iV + 7344] 2 



54(iV + 6) 3 VN 2 + 132N + 36 

The anomalous dimensions at these fix points give the critical exponents \ jv = 2 + 
Vm(9\ K) and j/v = 2- r] a (g*, XI) 



1 5N + 6 + y/N 2 + 132iV + 36 

v ~ ~ 6(7V + 6) C 



VN 2 + 132N + 36(3iV 3 + 109iV 2 + 510iV + 684) 

le 2 , (21) 



-3iV 4 - 658iV 3 - 333iV 2 - 15174iV + 4104] 2 



54(iV + 6)ViV 2 + 132iV + 36 
7_ N (UN + 6)VN 2 + 132N + 36 + 52iV 2 - 57iV + 36 2 

^- 2 "ivT6 e mWTW e ' 
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2.4 Gross-Neveu and Higgs- Yukawa Comparison 

The e = 4 — d expansions of the Gross-Neveu exponents || ||, [ToJ] are 

1, , 13 2 3 , ,1 

+ (39fc -^- »«C(|)-^ + ...)^, (23 ) 

- e + (fc -^-^ + ...)i 

51 „ 192C(3) + 281 , , 1 
+ (-36e+ ye 2 + ^^ e +'" ) ^' (24) 

while the 1/N expansion of Eqs. ([21] and ^2|) gives 

1 6 396 26136 

-lay = 2 - e +(-]y + ^-^r + --> 

, 13 1125 48951 . 2 



7, ,6 36 216 . 

, 7 51 1215 x 2 

Up to order e 2 and 1/iV 2 the two models agree as expected. In the GN 1/iV-expansion, 
the e 2 terms are comparable to the e ones and the e 3 is relatively small in 1/v and of the 
same magnitude in 'j/v. In contrast, the HY e-expansion shows that the coefficients 
of the 1/N expansion are always rapidly increasing, in particular in the case of 1/v. 
Thus a resummation for the GN v and the HY 7/1/ has to be made to improve the 
corresponding estimates. The necessity of such a resummation is also manifest from 
the importance, at low N, of the 1/Nj contribution for the GN u, and the e 2 one for 
the HY j/p (about 20% for JV=8). 

Because of the lack of information on asymptotic behaviour, we use a simple Pade- 



Borel resummation |14| with arbitrary choice of function, instead of the more sophis- 



ticated Borel resummation [|15|| . For an expansion 

A(x) = 1 + aix + a 2 x 2 + 0(x 3 ), (27) 

we write 

1 r°° 

A(x) = -J dte~ t/x [l - ait - (a 2 /2 - a\)t 2 Y . (28) 
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These formulae are directly used for the GN v obtained from Eq. (gj) [10], while we 
first expand the HY 7/z/ of Eq. (^) around the N = 00 point as 



j\ HY = (2 - e)(l + me + a 2 e 2 ) + 0(e 3 ), (29) 

and resum with Eq. ( p8|) only the second bracket. 

The comparison of the resulting critical exponents for the two 3d models as a function 
of the fermion number N is summarised in Fig. 1, where the dotted lines represent 
the computation of the GN \ jv and the HY 7/V without the resummation procedure 
described above. The difference between the two models is small except for v at low 
N. The data point at N = 48 comes from Ref. ||, while those at N = 8 result from 
the simulation described in the next section. 



3 Numercial Results 

Here we present our simulation and the analysis leading to estimates of the critical 
indices for the 3d Gross-Neveu model at iV = 8. 

3.1 Simulation of the lattice Gross-Neveu model 

For the numerical simulation we consider the effective action Eq. flrOp with Nl = 1 
which corresponds to N = 8, and use an exact Hybrid Monte Carlo algorithm. It has 
a point update of 8/xs-14/xs on a Cray Y-MP, increasing with lattice size. This is due 
to the fact that more conjugate gradient steps are needed to invert the fermion matrix 
for large lattices. We perform runs on lattice sizes 8 3 , 12 3 , 16 3 , 24 3 and 32 3 . We use 20 
time steps of length 0.2, except for the largest lattice where the time step is reduced to 
0.05. As a rule, measurements are carried out every 5th trajectory. Details regarding 
the runs are listed in Table |I[ The integrated autocorrelation time Tint quoted is that 
for (a 2 ). 

To analyse the data we use a variant of multihistogram reweighting analysis which 



does not require the binning of data [16]. This is used to obtain the values of observ- 
ables in between the simulated data points. These points are close enough, and the 
simulations long enough, to produce bosonic energy distributions that fill the whole 
coupling range of interest. 

3.2 The scaling and critical exponents from MC data 

Let us first consider the critical coupling and the critical exponent v, which describes 
the behaviour of the correlation length near the phase transition. We define the renor- 
malised coupling as 

9R = { -^y. (30) 
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Table 1: Statistics of the simulations. 



Size 




A 


Trajectories 


Tint 


8 3 





7875000 


110000 


8 


8 3 





8156250 


110000 


11 


8 


o 


8437500 


190000 


18 


12 3 





7875000 


90000 


8 


12 3 





8156250 


430000 


11 


12 3 





8437500 


90000 


15 


16 3 





7875000 


90000 


9 


16 3 





8156250 


282000 


11 


16 3 





8437500 


100000 


17 


24 3 





8156250 


432240 


14 


24 3 





8184375 


190000 


4 


32 3 





8167500 


267160 
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This expression for g R lacks a constant factor and is the inverse of the usual definition, 
but this does not affect its scaling properties. The scaling of qr is extremely simple: 

g R = f(L 1 H), (31) 

where / denotes a universal scaling function, L is the linear extent of the lattice and 
t = (A -1 — A c ~ )A C is the reduced coupling ("temperature"). The subscript c refers to 
the infinite volume critical coupling. 



We can determine A c by noting that, according to Eq. (|31~D, the curves of g R for 
different lattice sizes should cross at A c , up to scaling violations visible on too small 
lattices. The simulation results are shown in Fig. 2. The reweighting analysis gives 
crossings at A c = 0.820(2) for 8 3 and 12 3 , A c = 0.817(1) for 12 3 and 16 3 , A c = 0.815(1) 
for 16 3 and 24 3 and A c = 0.817(3) for 24 3 and 32 3 . We thus conclude that the GN 
model has a second order phase transition A c = 0.815(3). The renormalised coupling 
at A c is (g R ) c = 0.473(4). 

The usual way to avoid refering to A c in the critical exponents determination is to 
use thermodynamic quantities that peak in the scaling region. This is possible since, 
according to the scaling ansatz, L l l u t is constant at the maxima. Unfortunately, there 
are no quantities whose scaling behaviour provide a direct estimate of the exponent v. 
Therefore, if one tries to measure v one also has to specify the value of t, and thus A c . 

We can relax this requirement by noting that the scaling formula derived above is 
independent of the critical coupling and is valid in the whole critical region as long as 
the scaling violations can be neglected. Hence, we can perform a scaling analysis to 
quantities which do not necessarily peak at the critical coupling. Moreover, the latter 
need not even be specified: invering Eq. (0), L x l v t can be expressed as a function of 
gn and the finite size analysis can be made at constant value of Qr JT7| ]. We have to pay 
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a price, though, since measuring g R can be demanding. It is also crucial to eliminate 
t, since its uncertainty contributes a lot to the errors in the exponents. 

In order to extract the critical exponent u, we consider the logarithmic derivative of 
g R with respect to the reduced coupling 

D ee = L^Q^- ee L^F(L^t), (32) 

with F a new universal scaling function. Inverting the scaling equation of g R for L x l v t 
we obtain 

D = L^G(g R ), (33) 

where G is a scaling function. 

On the other hand, D is a correlator of powers of the a-field and the bosonic energy 
S = l/2J2n a n from the definition of average quantities : 

which can thus be measured numerically (i.e. by means of MC simulations). 
From this measurement, we determine v from a fit at constant g R to 

In D\ gR=const = \ jv In [L] + const (35) 

In Fig. 3, the value of v is shown as a function g R , together with the corresponding \ 2 
value of the fit. The scaling behaviour is realized for the entire fitting range 

X 2 < 0.5 (36) 

for three degrees of freedom. The errors on v are coming from a fit which uses the 
errors of the original data. These were obtained by a jackknifed reweighting analysis. 
The estimate obtained with g R = (g R ) c = 0.473 is 

v = 1.00(4). (37) 

Notice that the dependence on (g R ) is indead weak, and that we did not have to specify 
A c in our fits. To see how well our data is actually scaling we display g R versus tL l / u 
with our MC estimates of v and A c in Fig. 4. 

Using the hyperscaling relations, only two exponents are independent. As a second 
exponent we choose 2 — 77 = 7/V. This governs the behaviour of the susceptibility x 
near the critical point 

X=(a 2 )-(a) 2 = L^g(L^t), (38) 

where g denotes a scaling function. However, in numerical simulations there is a 
problem concerning the susceptibility: on finite lattices the average of the a field is 
always zero. 
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The use of absolute values in definition could distort the scaling behaviour and 
may lead to wrong exponents. To overcome this we used the susceptibility on the 
symmetric side [[Tj| A > A c where 

X=(a 2 )=L^g(L^t). (39) 
Eliminating again t as in the case for D (and g R ), we obtain 

In (& 2 )\ <m=const = - In L + const . (40 ) 
and get from a fit to the measured L dependence 

1 = 1.246(8). (41) 

Fig. 5 shows the results of the fit. Notice that the fit is valid for the whole critical 
range (x 2 < 0.4 for 3 d.o.f) and the deviation as a function of g R is very small as 
expected. The value we quote is taken at (<?ij) c = 0.473. From Fig. 6 one can see that 
the scaling of the data is excellent with our values of exponents: within error bars all 
the data from different lattice sizes lie on the same curve. 

As a check of consistency with hyperscaling, we can measure other critical exponents. 
The expectation value of the sigma field acts as an order parameter for the discrete 
chiral symmetry Eq. (Q), which is preserved on the lattice. On a finite lattice, the 
absolute value of a yields an estimate for the combination /3/v, which is shown in Fig. 
7. At g R = {g R ) c = 0.473 we get 

- = 0.877(4), (42) 
v 

with x 2 < 0-3 for 3 d.o.f. With this value of exponents the quality of scaling is again 
excellent, as one can see from Fig. 8. 

The combination /3/u is connected to j/u through the hyperscaling relation 

f = i(i-J). (43) 

Using the value of 7/1/ given in ([11]) this gives j3/u — 0.877(4), which is in complete 
agreement with estimate obtained above (f4~2"D. This agreement is noteworthy since we 
used the absolute value of a, which can lead to a distortion of the scaling relation. At 
least in our case we see that it does not. Also, the definition of the susceptibility with 
the absolute value of a leads to identical results. However, for the standard method of 
measuring the critical exponents from the scaling behaviour of thermodynamic quan- 
tities at their peak values, the usage of the absolute value (a) may result in a change 
in the position of the peaks and thus make the scaling analysis dubious. 

The heat capacity should give the combination a/ v. The hyperscaling relation pre- 
dicts a value of —1. This means that heat capacity does not diverge at the critical 
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Table 2: The critical exponents obtained from different methods. 
The numbers with a star are obtained with resummation. 



Exponent 


MC e e 2 


1/N f 


1/Nf 


V 


1.00(4) 0.9545 0.9480 


1.135 


0.903* 


l/v 


1.246(8) 1.4285 1.237* 


1.270 


1.2559 



point. In fact, it is dominated by its regular part which makes it impossible to extract 
a/ v. In order to do this we would need the second derivative of the heat capacity with 
respect to t. This quantity would diverge as t~~( a ' v+2 ' ~ t~ l . Unfortunately the quality 
of the MC data deteriorates as higher order derivatives of the free energy are taken: 
the 4th derivative is out of reach in the present simulation. 

All of the previous analysis relied heavily on reweighting the data from a finite set 
of couplings to a very dense set of couplings. This enabled us to accurately explore the 
dependence of the variables on the renormalised coupling g^. We note that both our 
method of analysis and the number of trajectories used allow us to achieve a better 
determination of the critical exponents than was achieved in a comparable analysis for 
N f = 12 §. 



4 Conclusions 

We have performed a high statistics simulation of the 3d GN model with two flavours of 
4-spinors. We show that it has a second order phase transition at A = 0.815(3). Hence, 
it leads to a continuum field theory, which is characterised by critical exponents which 
we have measured. This proves numerically that the GN model is renormalisable in 
three dimensions, even for a small number of flavours. The transition point is close 
to the 1/Nf expectation = 0.80, but the 1/Nj correction can be as significant as 
in 2d calculation jll| . Table |^ displays the results from our simulations together with 
estimates obtained by other methods: the e-expansion is for the HY model to first 
order by Zinn- Justin 0, and to second order as presented above, the 1/A^-expansion 
for the GN model calculated to one loop by Hands et al. || and to order 1/N 2 by 
Gracey M, [8], [10] and by Derkachov et al. ||. The second order contributions of the HY 



7/V and the GN v are large and the corresponding expressions have been resummed 
as explained in Sect. 2.4. 

The striking feature of the data is that the HY e-expansion results at the two-loop 
level are in very good agreement with the simulation values. Even without resumma- 
tion, which can be found quite arbitrary, the direct result of 7/z/ is not very far from 
the data point as seen in the Fig. 1. The GN second order l/iVy-expansion works very 
well for 7/z/ which has a small 1/Nj correction. Concerning the GN z/, even though 
the resummed value is not too far off from our numerical result, the discrepancy does 
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suggest that higher order terms may be important. However, the agreement with the 
HY v shows that no new phenomenon appears at small N. 

As a whole, our results strongly support the conjecture that these models are equiv- 
alent even in three dimensions, where they are not trivial, and that the properties 
inferred from perturbation theory are valid at low fermion number. 
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Figure captions 



Figure 1. Comparison of critical exponents obtained with different means as 
function of the effective fermion number N. Solid lines for the e 2 Higgs- Yukawa model, 
dashed lines for the 1/N 2 Gross-Neveu model, dotted lines for HY 7/V and GN v direct 
results (without resummation). The data point at N = 48 is for v from Ref. ||, those 
for N = 8 result from our simulation. 

Figure 2. The renormalised coupling g& as function of the coupling A for lattice 
sizes 8 3 , 12 3 , 16 3 , 24 3 and 32 3 in order of increasing slopes. The results of simulations 
without the reweighting are shown as circles. 

Figure 3. The critical exponent v as function of the value of g^. The corresponding 
X 2 plot gives the quality of the fit. 

Figure 4. The renormalised coupling g^ as function of tL l l u for different lattice 
sizes (8 3 is labeled with plus, 12 3 with octagons, 16 3 with squares, 24 3 with circles and 
32 3 with diamonds). The v and A c have the measured MC values. 

Figure 5. The critical exponent 7/V as function of the critical value of g^. The \ 2 
gives the quality of the fit. 

Figure 6. The combination a 2 L 3 ^^^ u and g^ as a function of tL 1 ^ for different 
lattice sizes (8 3 is labeled with plus, 12 3 with octagons, 16 3 with squares, 24 3 with 
circles and 32 3 with diamonds). The u, j/u and A c have the measured MC values. 

Figure 7. As in Fig. 5, but for the critical exponent /3/u. 

Figure 8. As in Fig. 6, but for the combination \a\LPI v . 
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